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NUMERICAL  PREDICTION  OF 


VORTEX  FLOWS 


F.  Boysan  and  J.  Swithenbank 

Department  of  Chemical  Engineering  and  Fuel  Technology 
University  of  Sheffield 


Abstract 

Numerical  predictions  of  turbulent  vortex  flow  in  cyclone  chambers  are 
made  with  the  aid  of  a  two-equation  model  of  turbulence  which  entails  the 
solution  of  elliptic  differential  equations  of  the  transport  of  turbulence 
energy  and  its  dissipation  rate.  The  n on-homogeneous  and  anisotropic  nature 
of  turbulence  is  accounted  for  by  employing  a  two-viscosity  concept. 

Reasonable  agreement  is  obtained  between  the  calculated  axial  and  tangential 
velocity  profiles  and  the  experimental  data  available  in  the  literature. 
Introduction 

Cyclone  dust  separators,  fluidic  vortex  amplifiers  and  swirl  combustors 
are  only  a  few  of  the  many  engineering  applications  of  the  cyclone  chamber. 

The  vortex  motion  in  this  device  is  created  by  the  tangential  introduction  of 
the  fluid  into  a  cylindrical  chamber  with  a  single  axial  exit  in  one  end  plate 
Numerous  configurations  are  possible  depending  on  the  number  and  position  of 
the  tangential  entry  ports,  the  location  of  the  axial  exit  and  the  length  of 
diameter  ratio  of  the  cylindrical  chamber.  Despite  its  simple  geometry,  the 
aerodynamics  of  the  cyclone  chamber  is  extremely  complicated.  A  typical  flow 
pattern  is  depicted  in  Fig.  1  which  is  characterized  by  several  annular  zones 
of  forward  and  reverse  flow,  entrainment  of  the  fluid  from  the  outside  at  the 
exit  and  the  high  degree  of  retention  of  swirl. 

In  the  past  there  have  been  several  experimental  studies  of  the  vortex 
structure  in  the  cyclone  chamber.  Baluev  and  Troyankin  [1]  have  reported 
measurements  of  the  all  three  components  of  velocity  in  23  designs  with  the 
aid  of  a  calibrated  five  channel  probe.  Hot-wire  anemometer  measurements  of 
not  only  the  mean  but  also  the  fluctuating  velocity  components  have  been  made 


by  Ustitnenko  and  Bukhman  [2j.  Observations  of  the  flow  pattern  in  a  cyclone 
chamber  using  smoke  have  been  conducted  by  Smith  [3],  who  also  measured  the 
axial  and  tangential  time  averaged  velocity  components  using  a  cobra  probe  to 
determine  direction  and  a  special  pilot-static  to  measure  the  magnitude. 
Unfortunately,  vortex  flows  are  extremely  sensitive  to  disturbances  created 
by  probes  [3,4,3]  and  hence  these  measurements  are  perhaps  of  questionable 
value.  Non-intrusive  optical  diagnostic  techniques  have  recently  started  to 
become  available  and  at  the  time  of  writing  of  this  paper  LDA  measurements  iu 
a  vortex  tube  have  been  reported  by  Escudier  et  al  [6], 

The  theoretical  analysis  of  swirling  flows  on  the  other  hand  is 
extremely  difficult  and  little  progress  in  this  area  has  so  far  been  achieved. 
Apart  from  the  fact  that  the  governing  equations  are  non-linear  and  strongly 
coupled,  the  structure  of  turbulence  is  non-homogeneous  and  anistropic, 
therefore,  not  suitable  for  analysis  by  closure  hypotheses  developed  for 
quasilinear  flows.  Numerical  solutions  of  laminar  flow  in  a  cyclone  chamber 
have  been  obtained  by  Lilley  and  Vatistas  [7].  However,  the  purpose  of  their 
efforts  have  been  to  demonstrate  the  capabilities  of  a  certain  numerical 
scheme  rather  than  to  arrive  at  quantitative  conclusions  about  the  structure 
of  the  flow. 

The  object  of  the  present  study  is  to  provide  numerical  predictions  of 
turbulent  vortex  flow  in  cyclone  chambers  and  to  investigate  the  applicability 
of  the  existing  models  of  turbulence. 

Mathematical  Formulation 

The  equations  required  for  the  description  of  the  aerodynamic  pattern  in 
cyclone  chambers  express  the  fluid  flow  balance  of  mass  and  momentum,  and  are 
given  here  in  cylindrical  co-ordinates  (x,  r,  6)  best  suited  to  the  geometry 
of  the  problem.  This  approach  is  supported  by  experimental  observations  that 
the  flow  in  cyclones  with  four  or  more  tangential  entries  loses  its  three 
dimensional  character  at  a  short  distance  from  the  ports  and  beeoaws  axially 
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symmetric.  The  derivatives  with  respect  to  the  0  co-ordinate  have  been 
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In  the  above  equations,  u,  v  and  w  are  the  components  of  velocity  in  the  x,  r 
and  0  directions  respectively,  p  is  the  pressure,  p  is  the  density  and  p  is 
the  viscosity. 

It  is  assumed  that  these  equations  are  valid  for  the  present  problem 
provided  that  all  the  flow  variables  and  fluid  properties  are  represented  by 
the  corresponding  time-mean  values,  p  is  now  the  effective  viscosity  which 
is  the  molecular  viscosity  augmented  by  the  turbulent  contribution.  The 
distribution  of  the  effective  viscosity  is  determined  by  the  ’turbulence 
model'  employed. 

Many  models  of  turbulence  exist,  some  involve  the  calculation  of  the 
effective  viscosity  from  a  prescribed  length  scale,  the  others  require  the 
solution  of  one  or  more  partial  differential  equations.  The  two-equation  k-e 
model,  which  is  of  moderate  complexity,  has  been  extensively  used  by  many 
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investigators  and  proven  to  be  adequate  over  a  wide  range  of  flow  situations. 
This  model  entails  the  solution  of  two  transport  equations  of  turbulence 
characteristics,  namely  the  local  energy  of  the  fluctuating  motion  K  and  the 
energy  dissipation  rate  e.  These  equations  can  be  written  as: 
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Where  a,,  and  a  are  the  turbulent  Prandtl  numbers  for  K  and  e,  and  G  is  the 
K  e 

rate  of  generation  of  K  which  is  defined  as: 


Solution  of  the  above  equations  for  K  and  e  allows  the  local  turbulence  and 
effective  viscosities  to  be  evaluated  from: 
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(8),  (9) 


respectively.  The  coefficients  C^,  C CQ,  C^,  oK  and  are  constants  which 
are  assigned  the  following  values,  =  1.44,  =  1.92,  Cp  =  1.0,  =  0.09, 

a  =  1.0  and  a  =  1.3  as  recommended  by  Launder,  Priddin  and  Sharraa  [8], 

Despite  the  success  with  which  the  K-e  model  given  by  equations  (5)  and 
(6)  has  allowed  prediction  of  numerous  turbulent  flows,  its  ability  to  repro¬ 
duce  the  strong  influence  of  swirl  has  not  been  sufficient.  It  has  been 
suggested  in  the  past  that  the  effects  of  rotation  can  be  accounted  for  by 
making  the  mixing  length  a  function  of  a  'swirl*  Richardson  number  [9,10]. 
This  basic  approach  is  not  compatible  with  the  K-e  model  where  the  distribu¬ 
tion  of  the  length  scale  is  obtained  from  the  solution  of  differential  trans¬ 
port  equations  rather  than  prescribed.  A  modification  consistent  with  the 
two  equation  model  has  been  proposed  by  Launder  et  al  [8],  which  involves  the 
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replacement  of  the  constant  in  the  transport  equation  of  the  energy 
dissipation  rate  by 

C2*  “  C2  (1  '  C3  Ri)  (10) 


Where  is  an  additional  constant  and  Ri  is  the  local  Richardson  number 
defined  as  follows 
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Furthermore,  in  view  of  the  recent  experimental  results  [ll]  which  dispute 
the  assumption  of  isotropy  on  which  the  K-e  model  rests,  a  simple  two- 
viscosity  approach  proposed  by  Lilley  [12]  is  adopted  in  the  present  study. 
The  effects  of  anisotropy  are  accomodated  in  the  conservation  of  momentum 
equation  in  the  0-direction  by  employing  an  effective  exchange  coefficient 
related  to  the  effective  viscosity  calculated  from  (9)  by 


u  cc  -  p  ,.,/a 
eff,w  eff  w 


where,  a  is  a  constant 
w 

Solution  Procedure 

The  equations  presented  in  the  preceding  section  are  elliptic  in  the  x 
and  r  directions  because  they  contain  the  second  order  derivatives  with  respect 
to  both  these  directions.  The  solution  procedure,  therefore,  needs  two 
dimensional  storage  and  requires  iteration.  Before  an  attempt  is  made  to 
solve  the  equations  however,  they  must  be  reduced  into  their  finite  difference 
analogues.  This  is  achieved  by  integrating  the  equations  over  a  computational 
cell.  The  resulting  algebraic  equations  can  be  represented  in  the  following 
common  form: 
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Where  the  A' 8  are  the  coefficients  which  contain  the  contributions  from  the 
convective  and  diffusive  fluxes,  and  are  the  components  of  the 
linearized  source  term.  The  set  of  algebraic  simultaneous  equations  (12)  are 
solved  by  the  SIMPLE  algorithm  of  Patankar  and  Spalding  [13],  This  algorithm 
involves  the  solution  of  the  momentum  equations  for  a  guessed  pressure 
distribution  to  give  a  first  estimate  of  the  velocity  fields.  Corrections 
to  the  pressu.e  and  velocity  fields  are  then  obtained  from  the  pressure 
equation  of  continuity.  These  corrections  are  such  that  the  resulting 
velocities  will  satisfy  continuity. 

Boundary  Conditions 

Because  of  the  elliptic  nature  of  the  equations,  a  complete  description 
of  the  problem  considered  necessitates  the  specification  of  the  boundary 
conditions  at  all  the  boundaries  of  the  domain  of  integration. 

In  the  near  wall  regions,  the  well  known  wall  functions  are  matched  with 
the  algebraic  equations  (12)  to  preclude  fine  grid  calculations  in  this 
region.  The  usual  practice  is  to  cut  the  link  between  the  boundary  and  near¬ 
wall  points  by  setting  the  appropriate  coefficient  to  zero,  and  to  insert  the 
wall  influence  by  way  of  the  linearized  source  terms.  The  specific  wall 
functions  employed  in  the  present  study  are 


T”  yp  x'1  In  (EYj) 


Y+  «  p  (K  C  *)*  -2- 

p  p  y  y 


*p  ‘  (Cn‘  V3/2/*yP 


Where,  t  is  the  shear  stress  at  the  wall,  u  ,  K  and  e  are  the  velocity, 
*w  P  P  P 

turbulence  energy  and  dissipation  rate  at  the  near  wall  node  respectively, 
y^  is  the  distance  from  the  wall,  x  and  E  are  constants. 
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The  boundary  conditions  at  the  axis  of  symmetry  are  of  the  zero  normal 
gradient  type  for  all  the  variables  except  the  radial  and  tangential 
components  of  velocity  which  are  themselves  zero  there. 

Although  the  conditions  at  the  inlet  are  specified  once  and  for  all  and 
do  not  need  updating  during  the  course  of  the  solution  procedure,  those  at 
the  exit  are  not  known  beforehand.  An  analytical  difficulty  associated  with 
confined  vortex  flows  is  their  extreme  sensitivity  to  downstream  boundary 
conditions  [ 14] .  In  order  to  mimimize  the  effects  of  the  downstream  condi¬ 
tions  on  the  flow  inside  the  cyclone  chamber,  the  domain  of  integration  is 
extended  further  downstream  of  the  chamber  exit.  The  boundary  conditions 
imposed  at  this  end  are  of  the  gradient  type  for  all  the  dependent  variables. 
Results 

The  calculations  were  performed  for  a  cyclone  chamber  of  diameter 
D  =  0.25  m,  length  L  =  0.386  m  and  exit  throat  diameter  Dc  =  0.1  m.  The 
fluid  was  assumed  to  be  fed  into  the  chamber  through  a  circumferential  slot, 
the  width  of  which  was  adjusted  so  that  the  flow  rate  and  the  inlet  velocities 
matched  the  experimental  conditions  reported  by  Ustimenko  and  Bukhman  [2]. 

The  ratio  of  tangential  to  radial  velocity  components  at  the  inlet  was  of  the 
order  of  10.  In  order  to  preserve  the  stability  of  the  numerical  solution 
procedure  and  to  obtain  convergence  for  such  high  degrees  of  swirl,  it  was 
found  necessary  to  reduce  the  under-relaxation  parameters  of  all  the  three 
velocity  components  to  0.25.  The  number  of  iterations  required  was  also  in 
direct  proportion  with  the  degree  of  swirl.  The  calculations  were  performed 
on  a  30  x  21  non-uniform  grid  in  the  x  and  r  directions  respectively. 

Fig.  2  shows  the  vector  plot  of  velocity  distribution  in  the  chamber 
predicted  using  the  conventional  K-e  model  of  turbulence.  It  can  be  seen 
th-t  the  incoming  flow  follows  the  chamber  wall  until  the  base  plate  is 
reached,  where  it  changes  direction  and  forms  a  reverse  stream.  The  fluid 
leaves  the  chamber  after  a  final  change  of  direction  without  sucking  any 


fluid  from  the  outside.  A  more  realistic  flow  pattern  is  shown  in  Fig.  3 
which  is  obtained  by  employing  the  two-viscosity  approach  including  the 
Richardson  number  correction  described  in  the  foregoing  sections.  It  is 
apparent  from  this  vector  plot  and  from  the  distribution  of  the  streamlines 
given  in  Fig.  4  that  the  incoming  stream  does  not  change  direction  in  the 
second  half  of  the  chamber,  but  the  streamlines  bend  radially  towards  the 
axis  of  symmetry  for  a  short  distance  before  they  become  almost  parallel  to 
it  and  flow  out  of  the  chamber.  This  time,  however,  fluid  is  sucked  into  the 
cyclone  chamber  from  the  outside,  which  is  in  accord  with  the  experimental 
observations  [1,2]. 

The  comparison  between  measured  (2)  and  calculated  axial  velocity  profiles 
are  displayed  in  Fig.  3.  it  is  interesting  to  note  that  all  the  calculated 
curves  which  correspond  to  the  K-e  model,  two-viscosity  approach  with  and 
without  Richardson  number  correction  show  fair  agreement  with  the  experimental 
data  in  the  region  0.5  K  <  r  <  R  ,  for  all  four  axial  locations.  In  the 
vicinity  of  the  axis,  however,  the  K-e  model  fails  to  predict  the  reverse 
flow  altogether.  The  two-viscosity  model  shows  some  reverse  flow  for  -  5 
and  =  0,  but  the  best  agreement  with  the  experimental  data  is  obtained  for 
a  =  2.5  and  =  10 

Fig.  6  shows  the  tangential  velocity  profiles  both  measured  and 
calculated  at  the  four  axial  locations.  It  is  apparent  that  the  K-e  model 
predictions  do  not  reproduce  even  the  qualitative  features  of  the  measured 
profiles  except  for  the  region  in  the  extreme  vicinity  of  the  cylindrical 
wall  and  underestimate  the  degree  of  retention  of  the  initial  swirl.  It  is 
found  that  better  qualitative  agreement  can  be  achieved  when  a is  increased. 
However,  as  seen  from  the  figure,  although  for  -  5  the  predicted  swirl 
velocities  are  well  above  the  measured  ones  there  is  not  enough  reverse  flow 


Finally,  Fig.  7  displays  the  profiles  of  the  kinetic  energy  of  turbulence. 

The  agreement  between  the  data  and  the  predicted  curves  is  indeed  poor.  It  is 

interesting  to  note  that  the  K-e  model  can  reproduce  the  experimental  profile 

for  r  >  0.5  R  ,  while  the  kinetic  energy  is  grossly  over-estimated  both  for 

-3 

0  *  5.0,  C0  =  0  and  a  =  2.5,  C.  =  10  .  This  suggests  that  a  also  may  be 

w  3  w  3  w  J 

a  function  of  Richardson  number  and  the  K-e  model  may  apply  after  all  in  the 
region  close  to  the  wall  of  the  chamber. 

Concluding  Remarks 

Numerical  predictions  have  been  presented  of  turbulent  vortex  flow  in 
cyclone  chambers.  The  turbulence  is  modelled  basically  by  the  usual  K-e  model. 
The  transport  equation  for  the  dissipation  rate  of  turbulence  energy,  however, 
was  modified  by  making  one  of  the  constants  a  function  of  the  local  Richardson 
number.  The  anisotropic  character  of  the  turbulence  was  accounted  for  by 
employing  a  different  effective  viscosity  for  the  0-direction  momentum 
equation.  Although  the  constants  a w  and  could  be  tuned  so  as  to  obtain 
reasonable  agreement  between  the  predicted  and  measured  axial  and  tangential 
velocity  profiles,  the  distribution  of  turbulence  energy  could  not  be 
reproduced. 

On  the  other  hand,  since  most  of  the  measurements  in  cyclone  chambers 
reported  in  the  literature  have  been  made  by  inseration  of  probes  into  the 
flow,  it  is  difficult  to  judge  whether  the  source  of  error  is  in  the  experi¬ 
ments  or  turbulence  modelling.  At  this  stage,  therefore,  one  must  conclude 
that  detailed  measurements  of  confined  vortex  flows  by  means  of  non-intrusive 
optical  techniques  are  required  urgently  to  validate  an  appropriate  turbulence 
model . 
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Figure  Captions 

Fig.  1  Flow  patterns  observed  in  a  cyclone  chamber  with  four 
tangential  entries  and  a  single  axial  exit  (from  [2]), 

Fig.  2  Velocity  vectors  in  the  cyclone  chamber  calculated  using 
the  k-e  model  of  turbulence. 

Fig.  3  Calculated  velocity  vectors  in  the  cyclone  chamber  for 

o  ■  2.5,  c-  =  10  3 
w  3 

Fig.  4  Calculated  distribution  of  streamlines  in  the  cyclone 
chamber  for  =  2<5  =  10“3 

w  3 

Fig.  5  Comparison  of  measured  and  predicted  axial  velocity 

profiles  at  four  downstream  stations.  Q  experiments, 

- k-e  model,  - k-e  model  with  0  =5  03=0. 0, 

- k-e  model  with  <^=2 . 5  C3=10“3 

Fig.  6  Comparison  of  measured  and  predicted  tangential  velocity 
profiles  in  the  cyclone  (for  legend  see  Fig.  5). 

Fig.  7  Comparison  of  measured  und  predicted  variation  of  kinetic 
energy  of  turbulence  (for  legend  see  Fig.  5). 


